
/**
* Use a NeutralDistribution to compute a NeutralDensity in the heliosphere
*
*  this is done by integrating over all velocity space for a desired grid
*
*  x, y, built around origin from + to - 'scale'.
*/
public class NeutralDensity {
	private static int gridsize = 40; // gridsize^2 spatial calcualtion points
	public static double AU = 1.49598* Math.pow(10,11); //meters
	private static double scale = 100; // edge size in AU
	public String filename = "testoutput3.dat";

	//private static double vMax = 10000000.0 // 10 k km/s


	// x and y are now the same but we hang onto them for clarity
	private double[] x = new double[gridsize];
	private double[] y = new double[gridsize];
	private double[] z = new double[gridsize*gridsize];

	private NeutralDistribution nd;
	private MultipleIntegration mi;


	/**
	*  Default, add lism params by hand here
	*/
	public NeutralDensity () {
		this(new NeutralDistribution(
			new KappaVLISM(
				new HelioVector(HelioVector.CARTESIAN,26000.0,0.0,0.0),
				100.0,8000.0,1.6),
			0.0, 0.0));//7*Math.pow(10.0,-8.0)));
	}


	/**
	*  Build this object from a distribution.
	*  0th order moment creation here
	*/
	public NeutralDensity(NeutralDistribution nd_) {
		nd = nd_;
		mi = new MultipleIntegration();
		//mi.setEpsilon(0.0001);
		nd.debug=false;

		//calculateGrid();
	}

	private void calculateGrid() {

		double delta = scale*AU/(double)gridsize;
		System.out.println("delta: " + delta);
		System.out.println("gridsize " + gridsize);
		int index = 0;
		for (int i=1; i<=gridsize; i++) {
			x[index]=0.0-(scale*AU/2)+delta*i;
			y[index]=0.0-(scale*AU/2)+delta*i;
			System.out.println("x["+index+"] : " + x[index]);
			index++;
		}

		index = 0;
		for (int i=0; i<x.length; i++) {
			for (int j=0; j<y.length; j++) {
				HelioVector pos = new HelioVector(HelioVector.CARTESIAN, x[i], y[j], 0.0);
				z[index]=vint(pos);
				index++;
				System.out.println("calculated z[" + (index-1) + "] : " + z[index-1]);
			}
		}

		file f = new file(filename);
		f.initWrite(false);
		for (int i=0; i<x.length; i++) {
			f.write(x[i]+"\n");
		}
		f.write("\n");
		for (int i=0; i<z.length; i++) {
			f.write(z[i]+"\n");
		}
		f.closeWrite();


		//JColorGraph jcg = new JColorGraph(x,y,z);
		//jcg.run();
	}


	/**
	* Do the integration!  Pass in a the position as a heliovector
	*/
	private double vint(HelioVector hp) {
		final HelioVector hpf = hp;
		FunctionIII dist = new FunctionIII () {
			// r,p,t here actually a velocity passed in...
			public double function3(double r, double p, double t) {
				return r*r*Math.sin(t)*nd.dist(hpf,new HelioVector(HelioVector.SPHERICAL, r, p, t));
			}
		};
		double[] limits = new double[6];
		limits[0]=0.0; limits[1]=100000;
		limits[2]=0.001; limits[3]=2*Math.PI;
		limits[4]=0.001; limits[5]=Math.PI;
		//limits[4]=Math.sqrt(2*nd.gmodMs/hp.getR()); limits[5]=300000;
		return mi.mcIntegrateSS(dist,limits, 100000);
	}

	/**
	* Use for testing first, then for starting routine to generate grid
	*/
	public static void main(String[] args) {
		// first we do some tests

		// first test - density at 1000AU,1000AU,1000AU

		/*HelioVector tv = new HelioVector(HelioVector.CARTESIAN, 100*AU, 100*AU, 100*AU);
		NeutralDistribution ndist = new NeutralDistribution(
						new GaussianVLISM(new HelioVector(HelioVector.CARTESIAN,50.0,0.0,0.0),
	        					  50000,8000.0)
	        		,0.0,0.0);

		NeutralDensity nd = new NeutralDensity(ndist);

	    System.out.println("test: " + nd.vint(tv));
		System.out.println("nd.counter: " + ndist.counter);  */

		NeutralDensity nd = new NeutralDensity();
		nd.calculateGrid();
	}
}